Cite this Article
Li Yuan-Lu, Meng Xiao, Ding Ya-Qing. Using wavelet multi-resolution nature to accelerate the identification of fractional order system. Chinese Physics B, 2017, 26(5): 050201  
Permissions
Using wavelet multi-resolution nature to accelerate the identification of fractional order system
Li Yuan-Lu1, 2, †, Meng Xiao2, Ding Ya-Qing2
B-DAT, School of Information and Control, Nanjing University of Information Science & Technology, Nanjing 210044, China
Jiangsu Collaborative Innovation Center on Atmospheric Environment and Equipment Technology, Nanjing University of Information Science & Technology, Nanjing 210044, China

 

† Corresponding author. E-mail: lyl_nuist@nuist.edu.cn

Abstract

Because of the fractional order derivatives, the identification of the fractional order system (FOS) is more complex than that of an integral order system (IOS). In order to avoid high time consumption in the system identification, the least-squares method is used to find other parameters by fixing the fractional derivative order. Hereafter, the optimal parameters of a system will be found by varying the derivative order in an interval. In addition, the operational matrix of the fractional order integration combined with the multi-resolution nature of a wavelet is used to accelerate the FOS identification, which is achieved by discarding wavelet coefficients of high-frequency components of input and output signals. In the end, the identifications of some known fractional order systems and an elastic torsion system are used to verify the proposed method.

PACS: 02.30.Jr;02.60.Cb;05.40.Ca;05.45.-a
Keyword:fractional order system;system identification;operational matrix;wavelet multi-resolution analysis
1. Introduction

In recent years, the FOS was widely used in various scientifc and engineering fields, such as physics,[1,2] chemistry,[3] biology,[4] signal processing,[5] control theory, etc.[6–8] The reason is that the fractional calculus can describe some phenomena much better than calculuses of their integer order counter-parts such as electrical-mechanical properties of materials,[9,10] viscoelastic material properties, etc.[11] However, how to identify an FOS from measured data of an actual system is still an open problem.

Currently, FOS identification methods mainly include time-domain identification methods and frequency-domain identification methods. In the time domain, the model parameters are estimated by minimizing the error between the output of the actual system and that of the identified system.[12–20] This method was first used to identify FOS by Poinot and Trigeassou[12] Subsequently, Jalloul et al.[13] used it to study the identification of the rotor skin effect in an induction machine. Malti et al. used the fractional Kautz-like basic function[14] and the optimal instrumental variable method[15] to identify the FOS. In addition, Liao et al.[16] and Wang et al.[17] used the subspace method to investigate commensurate FOS identification. In fact, the system identification can be expressed as an optimal problem, therefore, some intelligent optimization algorithms have also been utilized to identify FOS, such as particle swarm optimization algorithm,[18] genetic algorithm,[19] Luus–Jaakola algorithm, etc.[20]

In the frequency domain, Li and Yu[21] adopted the least squares algorithm to identify the non-integer order systems. Thomassin and Malti[22] developed a method to identify the state-space model of FOS, in which the parameters were determined by using the conventional subspace technique. Nazarian et al.[23] investigated the identification of FOS via frequency responses of the input and output of an actual system. Li et al.[24] combined the differential evolution algorithm and subspace identification algorithm to identify FOS with time delay. In Ref. [25], Levy’s identification method was extended to non-commensurable fractional system. In addition, orthogonal function is also used for FOS identification by Aoun et al.[26]

Recently, the operational matrix method based on orthogonal functions or polynomial series was diffusely used for analyzing the dynamic systems. The idea behind the method is that it converts the fractional calculus operator into a matrix. This method has been used to identify the FOS by Tang et al.[27] Li and Sun used the Block–Pulse operational matrix for analyzing the fractional order systems.[28] Other operational matrices include the Legendre operational matrix,[29–32] Chebyshev operational matrix,[33,34] Jacobi operational matrix,[35] operational matrices of triangular functions,[36] and Haar wavelet operational matrix.[37–39] However, for an actual system, the quantity of measurement data is large, which will lead to a high-dimensional operational matrix. As a result, the operational matrix method is inappropriate for identifying the fractional order system.

We observe that the wavelet coefficients of the input and output signal have a rapid downward trend. So, if you abandon certain high-frequency coefficients, you would find that the noise will be reduced and the dimensions of operational matrix will also be reduced, but the quality of the identification results does not significantly decrease.

2. Haar wavelet and operational matrix of fractional integration
2.1. Fractional calculus

The Riemann–Liouville fractional integration is defined as

where is the fractional integration order.

Compared with Riemann–Liouville definition, Caputo’s definition is convenient to describe the initial value problems of fractional differential equations, and it is defined as

where , , and γ is an integer.

Caputo’s integral definition has an important property

2.2. Haar wavelet operational matrix of fractional integration

The Haar wavelet is a kind of orthogonal function, which is defined as follows:

where

Any function can be expanded by Haar wavelets, i.e.,

where the Haar coefficients ri, , are determined by

In the practical application, the sum of infinite terms in Eq. (6) can be approximated with finite terms, namely,

where , the superscript T indicates the transposition, is the Haar coefficient vector, and is the Haar function vector.

We need N equations to obtain the coefficient vector . Hence, suitable collocation points are taken as

The N-square Haar matrix can be defined by

According to Ref. [28], we can know

where is the block pulse basis vector and is the block pulse operational matrix of the fractional integration, where block pulse functions may be defined in the interval as follows:
The in Eq. (11) is given by
where , , .

Let

where is called the operational matrix of the Haar wavelet of fractional integral of order α. From Eq. (10), we know that Haar matrix is an orthogonal function set, according to Ref. [40], the operational matrix of Haar wavelets can generally be calculated by using the block pulse function operational matrix

3. Multi-resolution analysis theory of the wavelet

Multi-resolution analysis means that the function is expressed as a combination of components with the time and frequency resolution. That is to say, a function is mapped to a series of nested approximation spaces.

A space, by two successive decompositions, can form a set of step by step contained subspaces. Spatial decomposition symbols are as follows:

where can be seen as the missed details that the low resolution function space ( is missed when they approach the high resolution function space ( , which means that is orthogonal complement space of in space . The sign “ ”denotes the orthogonal sum.

Each decomposition space has a set of orthogonal basis. Scaling function (approximation coefficients) is an orthogonal basis for , and Haar function (detail coefficients) is an orthogonal basis for . So, for any positive integer , space can be decomposed into

Thus, any function has the following multi-resolution representation:
Thus, can be divided into the following infinite orthogonal sum:
An arbitrary function can be developed as follows:
where , , and .

In practice, the infinite series in Eq. (20) is truncated, then it can be written as

where , , and .

According to Haar wavelet transform, we know that a signal can be approximated as

After wavelet transform, the coefficients of the signal are divided into the low-frequency part and the high-frequency part in each layer. Usually, the low-frequency part is regarded as the approximation of the original signal, especially in the case where the original signal is noisy. The low-frequency part has the effect of noise reduction. Thus, in each layer, the low-frequency coefficients are only half the coefficients of this layer. After the input and output signal are reconstructed through wavelet transform, the process of the system identification can be transformed into the identification of the wavelet coefficients. We discard the coefficients of high-frequency components layer by layer. Namely, we reduce the value of N proportionally, the quantity of data involved in system identification will be reduced by 2n, where n is the decomposition layer number. This paper is based on this idea to accelerate the fractional system identification.

In order to ensure that the signal is not distorted when we resample wavelet coefficients, the signal cannot be stratified at all times. That is to say, the value of N cannot keep down continuously. According to spectrum analysis of the original signal, we need to ensure that the new signal must satisfy sampling theorem when we abandon wavelet coefficients of high-frequency components of the signal. If the value of N is too small, the new signal will not represent the original signal accurately, the identification results will not make sense either. The final number of layers is determined by the error of the approximate components and the original signal. In general, the error is less than 95%.

Using different resolutions to represent the original input and output signal, a method for FOS identification is proposed. The advantage of the wavelet resolution feature is able to speed up the process of system identification.

4. Fractional system identification
4.1. Linear system identification

Considering a linear system described by the following differential equation:

where and are the input and the output of the system, respectively, bi is arbitrary positive number, and .

Let

By operating on both sides of Eq. (23) we can obtain the following integral equation:

Using Eq. (24) and operational matrix of integration, equation (25) can be rewritten as

Namely,

System identification is the process to obtain the identification parameters by the input and output of the system. Parameters , , , , , , and are unknown and to be identified, and these parameters are estimated by using the least-squares method. Firstly, we presuppose that the fractional derivative orders , , and are known, then we can obtain the model parameters , , , and .

Equation (27) can be written as

Let , , and .

Equation (27) can be rewritten as the following compact form:

The coefficient matrix can be obtained by the least-squares method

Then, varying fractional derivative orders , , and in a certain range, we can obtain the optimal derivative order by minimizing the error of Eq. (29) as the following:

4.2. Nonlinear system identification

Consider the following fractional order nonlinear system:

The derivation process is similar to that of the linear system. Assuming and , equation (32) can be simplified by a Haar wavelet operational matrix of integration as follows:

The Haar wavelet is a piecewise function and it can be extended with a block pulse function
Then
Let
Based on Eq. (12) and the definition of the block pulse function, we obtain that
Thus, combine Eqs. (36) and (37), equation (35) is simplified into

Like Eq. (38), we can obtain

According to Eq. (34), we have
So,
Substituting Eq. (41) into Eq. (33), equation (33) can be solved by
or

In system identification, and are known, namely, and in Eq. (43) are known. m, c, k, α 2, and α 1 are parameters to be identified.

Like linear system identification, the nonlinear system identification is presented briefly below.

Let

As a result, equation (43) may be rewritten in the following form:
Model coefficients can be solved by the least squares method

The following assessment criterion is used to find the optimal fractional derivative orders:

5. Illustrative examples

The effectiveness of the proposed method is illustrated by fractional order linear systems and nonlinear systems. The results of identification are assessed by the following formula:

where y is the output of the actual system and is the output of the identified system.

5.1. Identification of system under ideal case

Example 1 Consider a fractional order system, which is defined as

where , , , and is the unit step signal. Under zero initial condition, the output response of this system is obtained by the method introduced in subSection 4.1.

In the experiment, we set the change range and step size of the fractional derivative order for . The lengths of the input and output signal are both . Then we discard the wavelet coefficients of high-frequency components layer by layer. Namely, we take the lengths of wavelet coefficients to be , , , and , respectively, to identify the system. The identification results are shown in Table 1. The responses of the identified system for N = 64 and the actual system are shown in Fig. 1.

Fig. 1. (color online) Responses of the real model and the identified model of Example 1.
Table 1.

Identification results of Example 1 under different resolutions.

.

From Table 1, one can see that time needed for system identification is significantly reduced but the error is slightly increased.

Example 2 The second fractional order system[20] is described as

where , , , , and . This FOS was identified by a Luus–Jaakola algorithm in Ref. [20]. In order to compare our results with that estimated in Ref. [20], a random Binary signal is used as the input signal to excite the system. The input signal is shown in Fig. 2.

Fig. 2. (color online) Input signal of Example 2.

The ranges of fractional orders are set to be and . Identification results and the elapsed time of each experiment are recorded in Table 2. The identified results in Ref. [20] are also shown in Table 2. When , the output responses of the identified system and the real system are shown in Fig. 3.

Fig. 3. (color online) Responses of the real model and the identified model of Example 2.
Table 2.

Identification results of Example 2 under different resolutions.

.

From Table 2 and Fig. 3, one can find that when reducing the length of the wavelet coefficient, the running time of the computer is remarkably shorter but the real system can also be identified well.

Example 3 The third fractional order system is given by

We identify it as the following system . The identification results are listed in Table 3. The responses of the identified model ( ) and the real system are shown in Fig. 4.

Fig. 4. (color online) Responses of the real model and the identified model of Example 3.
Table 3.

Identification results of Example 3 identified as a second-order model.

.

This example indicates that the fractional order system can concisely represent a dynamic system.

Example 4 Consider a fractional order nonlinear system shown as follows:

where , , and . After the wavelet transform, we resample the wavelet coefficients by giving up the coefficients of high-frequency components to carry on the simulation. The lengths of the wavelet coefficients are , , , respectively. The identified results of the corresponding parameters and computer running time are shown in Table 4. The output responses of the real system and the identified system ( ) are shown in Fig. 5.

Fig. 5. (color online) Responses of the real model and the identified model of Example 4.
Table 4.

Identification results of Example 4 under different resolutions.

.

Example 5 A fractional doffing system is described by

where , , , , and , the fractional derivative orders vary in the following range: and . At first, the input and output signal are represented by term Haar wavelets and then we discard coefficients of high-frequency components. The identification results are given in Table 5. The output responses of the true system and the estimated system ( ) are shown in Fig. 6.

Fig. 6. (color online) Responses of the real model and the identified model of Example 5.
Table 5.

Identification results of Example 5 under different resolutions.

.

As can be seen from the result tables and response figures, whether fractional order linear system or nonlinear system is identified, the proposed method can achieve the desired identification parameters. Further, the identification time of each group listed in the tables indicates that using the multi-resolution property of the wavelet can greatly reduce the computer operating time.

The times of fractional linear and nonlinear system identification are shown in Figs. 7 and 8, respectively.

Fig. 7. (color online) Decomposition layer number and the consumed time.
Fig. 8. (color online) Decomposition layer number and the consumed time.

From Figs. 7 and 8, one can see that the consumption times of fractional linear and nonlinear system identification decrease as the increase of the decomposition layer. In Figs. 7 and 8, horizontal coordinates each denote the decomposition layer. When the decomposition layer is zero, the identification is directly carried out by using the coefficients of the input and the output of the system. When the decomposition layer is one, only the low-frequency coefficients are used to identify the FOS. Thus coefficients involved in identification are reduced by half of the original coefficients, so the consumed time for system identification is reduced.

5.2. Identification of systems under noise case

Wavelet multi-resolution characteristics can not only speed up the process of the system identification, but also have the ability to resist interference. The above five experiments about parameter identifications of FOS are carried out in the absence of interference. However, the actual signal often contains noise. In order to verify that the mentioned approach is immune to noise, we discuss the system identification in the case of noise. Considering Example 1, we add white Gaussian noise to output signals, and the signal-noise-ratios (SNRs) are 20 dB, 30 dB, and 40 dB, respectively.

In these experiments, the output responses with noise have been pre-processed by the wavelet de-noising method.

Under SNR = 20 dB, 30 dB, 40 dB, respectively the identification results are listed in Tables 6–8. Figure 9 shows the responses of the real system and the identified system with the largest identification error (N = 1024) when SNR is 30 dB.

Fig. 9. (color online) Under SNR = 30 dB, responses of the real model and the identified model of Example 1.
Table 6.

Under SNR = 20 dB, identification results of Example 1.

.
Table 7.

Under SNR = 30 dB, identification results of Example 1.

.
Table 8.

Under SNR = 40 dB, identification results of Example 1.

.
5.3. Identification of a multi-mass elastic torsion system

The following example is about the identification of a multi-mass elastic torsion system. Two cylindrical masses are connected to both ends of an elastic spring, and their central axis is fixed at the same level. A DC motor is connected to the left of the first mass. The DC motor drives the first mass, through the spring, which can make the second mass rotate. The encoder is used to collect the rotational speed of the second mass.

Input and output signals of the real system are shown in Fig. 10.

Fig. 10. (color online) (a) Input and (b) output signals of the multi-mass elastic torsion system.

Assume the identified model as follows:

Using the proposed method, the identified parameters of the FOS model and the IOS model are listed in Table 9. The error between the output of the real system and that of the identified system is also given in Table 9. The identified results are shown in Figs. 11 and 12.

Fig. 11. (color online) Multi-mass elastic torsion system identified as an FOS.
Fig. 12. (color online) Multi-mass elastic torsion system identified as an IOS.
Table 9.

Identification results of the multi-mass elastic torsion system.

.

From Figs. 11 and 12, we can see that the IOS and FOS both can model the multi-mass elastic torsion system, but the fractional order model is more accurate than the integral model.

6. Conclusions

A fractional order system identification method based on wavelet multi-resolution analysis is proposed. The method reduces the lengths of the input signal and the output signal by resampling wavelets coefficients of the input signal and the output signal, which can improve the efficiency of fractional order system identification. Numerical simulations, including some known fractional order system and the multi-mass elastic torsion system, confirm the validity of the above method.

Reference
[1] Zhang R X Yang S P 2011 Chin. Phys. 20 090512
[2] Zhang D X Zhao D Guan Z H Wu Y H Chi M 2016 Phys. A: Stat. Mech. Appl. 461 299
[3] Oldham K B 2010 Adv. Eng. Softw. 41 9
[4] Magin R L 2010 Comput. Math. Appl. 59 1586
[5] Li Y Pan C Meng X Ding Y Chen H 2015 Meas. Sci. Rev. 15 101 http://dx.doi.org/10.1515/msr-2015-0015
[6] Li T Z Wang Y Luo M K 2014 Chin. Phys. 23 080501
[7] Chen J Guan Z H Li T Zhang D X Ge M F Zheng D F 2015 Neurocomputing 168 698
[8] Chen J Guan Z H Yang C Li T He D X Zhang X H 2016 J. Frankl. Inst. 353 1672
[9] Cao H L Deng Z H Li X Yang J Qin Y 2010 Int. J. Hydrogen Energy. 35 1749
[10] Chen L F Cheng L Y 2014 Sci. Technol. Eng. 14 54 http://dx.doi.org/10.3969/j.issn.1671-1815.2014.32.011
[11] Yao G J Lv W G Song R L Cui Z W Zhang X L Wang K X 2010 Chin. Phys. 19 074301
[12] Poinot T Trigeassou J C 2004 Nonlinear Dyn. 38 133
[13] Jalloul A Trigeassou J C Jelassi K Melchior P 2011 Int. J. Comput. Sci. Iss. 8 801 http://dx.doi.org/10.1007/s11071-013-0833-8
[14] Malti R Aoun M Oustaloup A 2004 Proceedings of the 30th Chinese Control Conference March 21–24 2004 Hammamet, Tunisia 835 839
[15] Malti R Victor S Oustaloup A 2008 17th IFAC World Congress July 2008 Seoul, Korea 14379 14384
[16] Liao Z Peng C Wang Y 2011 30th Chinese Control Conference July 22–24 2011 Yantai, China 1636 1640
[17] Wang Y Liao Z Peng C Liang S Zhu Z T 2013 Control. Decis. 28 67 http://www.kzyjc.net:8080/CN/abstract/abstract12205.shtml
[18] Li X Han Z 2011 Control. Decis. 26 1627 http://www.kzyjc.net:8080/CN/abstract/abstract10357.shtml
[19] Zhou S X Cao J Y Chen Y Q 2013 Entropy 15 1624
[20] Li D Z Yu Z X 2008 J. Tsinghua University (Sci & Tech) 48 1742 (in Chinese) http://dx.doi.org/10.3321/j.issn:1000-0054.2008.10.010
[21] Li Y L Yu S L 2007 Acta Autom. Sin. 33 882 http://dx.doi.org/10.1360/aas-007-0882
[22] Thomassin M Malti R 2009 15th IFAC Symposium on System Identification July 2009 Saint Malo, France TuA3 6
[23] Nazarian P Haeri M Tavazoei M S 2010 ISA Trans. 49 207
[24] Li W Peng C Wang Y 2011 Int. J. Control Autom. Syst. 9 310
[25] Valério D Tejado I 2015 Signal Process. 107 254
[26] Aoun M Malti R Levron F Oustaloup A 2007 Automatica 43 1640
[27] Tang Y G Liu H F Wang W W Lian Q S Guan X P 2015 Signal Process. 107 272
[28] Li Y L Sun N 2011 Comput. Math. Appl. 62 1046
[29] Jafar H Yousefi S A Firoozjaee M A Momani S Khalique C M 2011 Comput. Math. Appl. 62 1038
[30] Bhrawy A H Doha E H Ezz-Eldien S S Abdelkawy M A 2016 Calcolo 53 2 http://dx.doi.org/10.1007/s10092-014-0132-x
[31] Bhrawy A H Alghamdi M A 2013 Adv. Differ. Equations 1 2 http://dx.doi.org/10.1186/1687-1847-2013-307
[32] Bhrawy A H Taha T M Alzahrani E O Baleanu D Alzahrani A A 2015 PloS One 10 e0126620
[33] Li Y L 2010 Commun. Nonlinear Sci. Numer. Simul. 15 2284
[34] Seifollahi M Shamloo A S 2013 Word Applied Program 3 85 https://www.researchgate.net/publication/262377450
[35] Doha E H Bhrawy A H Ezz-Eldien S S 2012 Appl. Math. Model. 36 4931
[36] Asgari M 2015 Iaeng Int. J. Appl. Math. 45 85 https://www.researchgate.net/publication/281951567
[37] Li Y L Zhao W W 2010 Appl. Math. Comput. 216 2276 http://dx.doi.org/10.1016/j.amc.2010.03.063
[38] Gao G Gu S 2008 J. Tsinghua University (Sci & Tech) 48 1821 (in Chinese) http://dx.doi.org/10.3321/j.issn:1000-0054.2008.10.026
[39] Li Y L Meng X Zheng B C Ding Y Q 2015 ISA Trans. 59 79
[40] Wu J L Chen C H Chen C F 2001 IEEE Trans. Circ. Syst I: Fund. Theor. Appl. 48 120